
## 2.5作图合集
# 加载所需库------
if (!require("pacman")) install.packages("pacman")
pkgs <- c("dplyr", "ggplot2", "readxl", "tidyr", "lubridate", "janitor", "stringr", "writexl",
          "psych", "magrittr", "sf", "spdep", "tmap", "ggspatial","RColorBrewer","here","scales")
pacman::p_load(pkgs, character.only = TRUE)

source("./code/2.3 人口数估发病率.R")


### 主图1：2005-2023年逐年逐月发病数------
yearly_month_incidence <- rate_month %>% 
  dplyr::select(region,year,month,incidence_rate = rate,cases = count, pop = value)

# 生成完整日期序列并填充缺失值
yearly_month_incidence1 <- yearly_month_incidence %>%
  mutate(date = ymd(paste(year, month, "01", sep = "-"))) %>%
  complete(date = seq(ymd("2005-01-01"), ymd("2023-12-01"), by = "1 month")) %>%
  mutate(cases = coalesce(cases, 0L))


# 步骤1：数据预处理
province_data <- yearly_month_incidence1 %>% 
  filter(region == "province") %>%  # 筛选省级数据
  mutate(date = ymd(date)) %>% 
  arrange(date)

# 计算双轴比例系数（关键）
scaling_factor <- max(province_data$cases, na.rm = TRUE) / max(province_data$incidence_rate, na.rm = TRUE)

# 步骤2：构建双轴图
ggplot(province_data, aes(x = date)) +
  # 主y轴：发病数柱状图
  geom_col(aes(y = cases, fill = "Cases"), 
           width = 20, 
           alpha = 0.7,
           linewidth = 0.3) +
  # 次y轴：发病率折线图
  geom_line(aes(y = incidence_rate * scaling_factor, color = "Incidence Rate"), 
            linewidth = 0.8,
            key_glyph = "path") +
  # 双轴标度设置
  scale_y_continuous(
    name = "Number of Cases",
    labels = label_number(scale_cut = cut_short_scale()),
    sec.axis = sec_axis(
      ~./scaling_factor,
      name = "Incidence Rate (per 100,000)",
      labels = number_format(accuracy = 0.1)
    ),
    expand = expansion(mult = c(0, 0.05))
  ) +
  # X轴时间刻度
  scale_x_date(
    breaks = seq(ymd("2005-01-01"), ymd("2023-01-01"), by = "2 years"),
    date_labels = "%Y",
    expand = expansion(add = c(30, 60))
  ) +
  # 颜色标度（符合IDP期刊风格）
  scale_fill_manual(
    name = "",
    values = c("Cases" = "#2B8CBE"),
    labels = c("Cases"),
    guide = guide_legend(order = 1)
  ) +
  scale_color_manual(
    name = "",
    values = c("Incidence Rate" = "#D62728"),
    labels = c("Incidence Rate"),
    guide = guide_legend(order = 2)
  ) +
  # 期刊风格主题
  theme_minimal() +
  theme(
    text = element_text(family = "Arial", color = "#333333"),
    panel.grid.major = element_line(color = "grey90", linewidth = 0.25),
    panel.grid.minor = element_blank(),
    axis.line = element_line(color = "black", linewidth = 0.4),
    axis.title.y.left = element_text(size = 10, face = "bold", margin = margin(r = 10)),
    axis.title.y.right = element_text(size = 10, face = "bold", margin = margin(l = 10)),
    axis.text = element_text(size = 8),
    legend.position = "top",
    legend.spacing.x = unit(0.3, "cm"),
    legend.key.width = unit(1.5, "cm"),
    legend.text = element_text(size = 9),
    plot.margin = margin(15, 20, 10, 15)
  ) +
  labs(x = NULL)

# 保存高清图
ggsave("./主图/Fig.1 2005-2023年逐年逐月发病数.png",dpi = 600)

### 主图1.2:年度累计热力图：月发病数-----


yearly_incidence <- readxl::read_excel( "./output/表格2.1_发病数和发病率_年度.xlsx")


### 绘图：整个新疆地区2005-2023年年度发病数和发病率（双Y轴图）

# 绘制双Y轴图
ggplot(yearly_incidence) +
  geom_col(aes(x = year, y = cases), fill = "#8DD3C7", alpha = 0.6) +
  geom_line(aes(x = year, y = incidence_rate * 1000 / 4), color = "#FB8072", size = 1) +
  geom_point(aes(x = year, y = incidence_rate * 1000 / 4), color = "#FB8072") +
  scale_y_continuous(
    name = "Case Count",
    sec.axis = sec_axis(~ . / 1000 * 4, name = "Incidence Rate (per 100,000)")
  ) +
  labs(title = "Annual Case Count and Incidence Rate in Xinjiang (2005-2023)", x = "Year") +
  theme_minimal() +
  theme(
    axis.title.y = element_text(color = "#8DD3C7"),
    axis.title.y.right = element_text(color = "#FB8072"),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )
ggsave("./主图/图2.2 整个新疆地区2005-2023年年度发病数和发病率.png")

### 主图1.3：月度累计热力图：月发病数------


monthly_cases_ratio <- readxl::read_excel("./output/表格2.1_发病数和发病率_月份累积.xlsx")

# 绘制双Y轴图
# 计算两个轴的最大值，确保次Y轴能够正确显示
max_cases <- max(monthly_cases_ratio$total_monthly_cases)
max_rate <- max(monthly_cases_ratio$incidence_rate)
scale_factor <- max_cases / max_rate

ggplot(monthly_cases_ratio) +
  geom_col(aes(x = month, y = total_monthly_cases), fill = "#8DD3C7", alpha = 0.6) +
  geom_point(aes(x = month, y = incidence_rate * scale_factor), color = "#FB8072") +
  geom_line(aes(x = month, y = incidence_rate * scale_factor), color = "#FB8072", size = 1,group = 1) +
  scale_y_continuous(
    name = "Case Count",
    sec.axis = sec_axis(~ . / scale_factor, name = "Incidence Rate (per 100,000)")
  ) +
  labs(title = "Monthly Case Count and Incidence Rate in Xinjiang (2005-2023)", x = "Month") +
  theme_minimal() +
  theme(
    axis.title.y = element_text(color = "#8DD3C7"),
    axis.title.y.right = element_text(color = "#FB8072"),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )


# 保存图表
ggsave("./主图/图2.3 整个新疆地区2005-2023年各月发病数和发病率.png")



# 全局参数设置 (核心配置)
color_breaks <- c(0, 500, 1000, 1500, 2000, 2500)  # 强制固定色阶
color_palette <- colorRampPalette(c("#FEE5D9", "#FB6A4A", "#CB181D"))(length(color_breaks)-1)
color_labels <- paste0(color_breaks[-length(color_breaks)], "-", color_breaks[-1])


# 函数：标准化热力图模板

idp_heatmap <- function(data, x_var, title) {
  ggplot(data, aes(x = factor({{x_var}}), y = nameE, 
                   fill = cut(total_cases, 
                              breaks = color_breaks,
                              include.lowest = TRUE,
                              right = FALSE))) +
    geom_tile(color = "white", linewidth = 0.3) +
    scale_fill_manual(
      values = color_palette,
      labels = color_labels,
      name = "Total Cases",
      drop = FALSE  # 强制显示所有色阶
    ) +
    labs(x = NULL, y = "City", title = title) +
    theme(
      # 基础设置
      text = element_text(family = "Arial"),
      panel.background = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      # 坐标轴
      axis.text.x = element_text(
        angle = 90,
        vjust = 0.5,
        hjust = 1,
        size = 10,
        color = "black",
        face = "bold"
      ),
      axis.text.y = element_text(
        size = 10,
        color = "black",
        face = "bold"
      ),
      axis.title.y = element_text(
        size = 12,
        margin = margin(r = 10),
        face = "bold"
      ),
      # 标题
      plot.title = element_text(
        size = 14,
        hjust = 0.5,
        margin = margin(b = 15),
        face = "bold"
      ),
      # 图例
      legend.position = "bottom",
      legend.title = element_text(size = 10, face = "bold"),
      legend.text = element_text(size = 8),
      legend.key.size = unit(0.5, "cm"),
      # 网格线
      panel.grid = element_blank()
    ) +
    coord_fixed(ratio = 0.8)
}

# 主图2.1：月度热力图
monthly_cases_by_city <- read_excel("./output/表格2.4.1_发病数和发病率_城市.xlsx") %>%
  mutate(month = factor(month, levels = 1:12, labels = month.abb));14*12  # 转换为月份缩写

(idp_heatmap(monthly_cases_by_city, month, "Monthly Brucellosis Cases (2005-2023)") +
    scale_x_discrete(breaks = month.abb))

ggsave("./主图/Fig2.1_Monthly_Heatmap.png", width = 8, height = 6, dpi = 300)

# 主图2.2：年度热力图

yearly_cases_by_city <- read_excel( "./output/表格2.4.1_发病数和发病率_各城市各月份.xlsx");
14*19

(idp_heatmap(yearly_cases_by_city, year, "Annual Brucellosis Cases (2005-2023)"))

ggsave("./主图/Fig2.2_Annual_Heatmap.png", width = 8, height = 6, dpi = 300)

#### 主图3：分阶段发病数和发病率--------

# 读取地图
# nameCE_county$county_code  <-  as.character(nameCE_county$county_code)
map_data <- st_read("./data/origin/新疆维吾尔自治区.json") %>%
  st_transform(crs = 4326) %>%
  mutate(NAME = name, GB1999 = as.character(adcode), city_code = substr(GB1999,1,4)) 

# 检查哪些几何对象是无效的
invalid_index <- which(!st_is_valid(st_geometry(map_data)))

if (length(invalid_index) > 0) {
  cat("发现无效几何对象，尝试修复...\n")
  
  # 对无效几何对象进行修复
  map_data$geometry[invalid_index] <- st_make_valid(st_geometry(map_data)[invalid_index])
  
  # 再次检查是否仍然有无效几何对象
  invalid_after_first_fix <- which(!st_is_valid(st_geometry(map_data)))
  
  if (length(invalid_after_first_fix) > 0) {
    cat("仍有无效几何对象，尝试使用缓冲处理修复...\n")
    
    # 使用 st_buffer 修复剩余的无效几何对象
    buffer_dist <- 0.000001  # 使用一个非常小的正值避免重复边问题
    map_data$geometry[invalid_after_first_fix] <- st_buffer(
      st_geometry(map_data)[invalid_after_first_fix],
      dist = buffer_dist
    )
    
    # 再次检查是否仍然有无效几何对象
    invalid_after_buffer <- which(!st_is_valid(st_geometry(map_data)))
    
    if (length(invalid_after_buffer) > 0) {
      cat("仍有无效几何对象，尝试简化几何对象...\n")
      
      # 尝试简化几何对象
      simplify_tolerance <- 0.0001  # 根据需要调整简化容差
      map_data$geometry[invalid_after_buffer] <- st_simplify(
        st_geometry(map_data)[invalid_after_buffer],
        dTolerance = simplify_tolerance
      )
      
      # 最后再检查一次有效性
      final_invalid <- which(!st_is_valid(st_geometry(map_data)))
      if (length(final_invalid) > 0) {
        warning("仍有无效几何对象无法自动修复，请手动检查这些对象。\n")
      } else {
        cat("所有几何对象已成功修复。\n")
      }
    } else {
      cat("所有几何对象已通过缓冲处理成功修复。\n")
    }
  } else {
    cat("所有几何对象已通过 st_make_valid 成功修复。\n")
  }
} else {
  cat("所有几何对象都是有效的，无需修复。\n")
}

# 确保 map_data 的几何列是最新的
map_data <- st_as_sf(map_data)

# 检查哪些几何对象是无效的
invalid_index <- which(!st_is_valid(st_geometry(map_data)))

# 打印结果以确认
print(head(map_data))

merged_data <- readRDS("./data/processed/merged_data.rds")


# 3 主图3：分阶段发病率地图--------

time_periods <- list(
  y2005_2023 = 2005:2023,
  y2005_2009 = 2005:2009,
  y2010_2015 = 2010:2015,
  y2016_2019 = 2016:2019,
  y2020_2023 = 2020:2023
)


# 修正后的发病率计算函数
calculate_true_incidence <- function(data, years) {
  data %>%
    filter(year %in% years) %>%
    group_by(county_code, NAME,nameE) %>%
    summarise(
      total_cases = sum(count, na.rm = TRUE),      # 时段总发病数
      total_pop = sum(pop, na.rm = TRUE),          # 时段总人口数
      avg_incidence_rate = (total_cases / total_pop) * 100000,  # 每10万人发病率
      .groups = 'drop'
    ) %>%
    mutate(avg_incidence_rate = ifelse(is.nan(avg_incidence_rate), 0, avg_incidence_rate)) %>%
    select(-total_cases, -total_pop)  # 可选：移除中间计算列
}

# 计算各时段真实发病率
incidence_rates <- map(time_periods, ~calculate_true_incidence(merged_data, .x))

# 命名结果列表
names(incidence_rates) <- names(time_periods)

# 将 incidence_rates 列表中的每个数据框合并到一个单独的数据框中，并且为每个数据框添加一个标识时间段的列
# 为每个子数据框添加时间段标识，并合并所有子数据框
combined_data <- map_df(incidence_rates, ~ .x %>% 
                          mutate(time_period = cur_group()), .id = "time_period")%>%
  mutate(time_period = factor(time_period, levels = names(incidence_rates)),
         incidence_rate = ifelse(is.na(avg_incidence_rate), 0, avg_incidence_rate))

writexl::write_xlsx(combined_data, "./output/发病率_分阶段.xlsx");5*(96)

# 查看合并后的数据框前几行
print(head(combined_data))

summary(combined_data$incidence_rate)

# 定义全局的 color_scale 范围
global_color_scale <- seq(0, 300, by = 50)

# 定义色板
palette <- c(brewer.pal(9, "OrRd")[c(1:7)]) #display.brewer.pal(9, "OrRd")
# palette <- brewer.pal(7, "Spectral") %>% rev() #display.brewer.pal(7, "Spectral")

# palette <- brewer.pal(7, "Spectral") %>% rev() #display.brewer.pal(9, "RdYlBu") #dev.new()


# 绘制地图的函数 
plot_incidence_map_tmap <- function(rate_data, title, global_color_scale, palette) {
  # 合并地图数据与发病率数据
  rate_map_data <- map_data %>%
    left_join(rate_data, by = c("GB1999" = "county_code"))
  
  # 创建 tmap 地图
  tm_shape(rate_map_data) +
    tm_polygons(col = "avg_incidence_rate",
                style = "fixed",
                palette = palette,
                breaks = global_color_scale,
                title = "Incidence (1/100,000)",
                # na.color = "white",
                textNA = "Data Missing",  # 图例中显示缺失值标签
                na.color = "white") +# 缺失值区域颜色设为白色
    tm_layout(
      main.title = title,
      main.title.position = "center",
      legend.position = c("right", "bottom"),
      legend.text.size = 0.8,
      legend.title.size = 0.9,
      frame = FALSE
    ) +
    tm_compass(position = c("left", "top"), size = 2) +  # 添加指北针
    tm_scale_bar()        # 添加比例尺
}

# 确保目录存在，如果不存在则创建
# dir.create("figs/逐年发病地图", showWarnings = FALSE, recursive = TRUE)

# sub("^y", "", "y2005_2009") %>% 
#   strsplit("_") %>% unlist() %>% as.numeric() %>% .[1]

# 分别绘制不同时间段的地图并保存图片
plots <- lapply(names(incidence_rates), function(period) {
  p <- plot_incidence_map_tmap(
    incidence_rates[[period]],
    title = paste("Avaraged Annual Incidence Rate for Peroid ", 
                  sub("^y", "", period) %>% 
                    strsplit("_") %>% unlist() %>% as.numeric() %>% .[1],
                  "-",
                  sub("^y", "", period) %>% 
                    strsplit("_") %>% unlist() %>% as.numeric() %>% .[2]),
    global_color_scale = global_color_scale,
    palette = palette)
  filename <- file.path("./主图", paste0("incidence_rate_", period, ".png"))
  tmap_save(p, filename = filename, dpi = 300)})


# 4 主图4：时空扫描地图--------
st_clean <- readxl::read_excel("./data/origin/st_result/st_clean.xlsx")

county_bt <- nameCE$county_code[nameCE$region %in% "county兵团"]

# 定义绘图函数
plot_cluster_map <- function(cluster_data, 
                             map_data, 
                             target_year,
                             title_prefix = " ") {
  
  # 过滤当年数据
  annual_cluster <- cluster_data %>%
    filter(year == target_year) %>%
    mutate(county_code = as.character(county_code))
  # map_data$adcode %>% unique()
  # 合并地理数据并处理缺失值
  cluster_map_data <- map_data %>% filter(!adcode %in% county_bt) %>%
    mutate(GB1999 = as.character(GB1999)) %>%
    left_join(annual_cluster, by = c("GB1999" = "county_code")) %>%
    mutate(Cluster_type = case_when(
      is.na(Cluster_type) ~ "Non-cluster",
      TRUE ~ Cluster_type
    )) %>%
    mutate(Cluster_type = factor(Cluster_type,
                                 levels = c("Primary Cluster", "Secondary Cluster", "Non-cluster")))
  
  # 期刊推荐配色（ColorBrewer 8-class RdYlBu反向）
  cluster_colors <- c(
    "Primary Cluster" = "#d73027",  # 深红色
    "Secondary Cluster" = "#fc8d59", # 橙色
    "Non-cluster" = "white"
  )
  
  # 创建地图
  base_map <- tm_shape(cluster_map_data) +
    tm_polygons(
      col = "Cluster_type",
      palette = cluster_colors,
      title = "Cluster Category",
      style = "cat",
      border.col = "gray40",
      lwd = 0.4,
      alpha = 0.9,
      legend.show = TRUE
    ) +
    tm_layout(
      main.title = paste0(title_prefix, target_year),
      # main.title.size = 1.2,
      main.title.position = "center",
      main.title.fontface = "bold",
      fontfamily = "Arial",
      legend.position = c("right", "bottom"),
      legend.bg.color = "white",
      legend.bg.alpha = 0.8,
      legend.frame = FALSE,
      # legend.text.size = 0.85,
      # legend.title.size = 1.0,
      frame = FALSE,
      outer.margins = c(0.02, 0, 0.02, 0)
    ) #+
    tm_compass(position = c("left", "top"), size = 2) +  # 添加指北针
    tm_scale_bar()
  
  return(base_map)
}

# 主程序 
# 设置输出目录
output_dir <- "主图/时空扫描地图"
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)

# 预处理数据（英文转换）
st_clean_en <- st_clean %>%
  mutate(
    Cluster_type = case_when(
      Cluster_type == "cluster I" ~ "Primary Cluster",
      Cluster_type == "cluster II" ~ "Secondary Cluster",
      TRUE ~ "Non-cluster"
    )
  )

# 生成年份序列（2005-2023）
years <- st_clean_en %>%
  pull(year) %>%
  unique() %>%
  sort()

# 批量导出地图
walk(years, ~{
  cat("Generating map for year", .x, "...\n")
  
  cluster_map <- plot_cluster_map(
    cluster_data = st_clean_en,
    map_data = map_data,
    target_year = .x
  )
  
  output_file <- file.path(
    output_dir,
    paste0("ClusterMap_", .x, ".tiff")
  )
  
  tmap_save(
    tm = cluster_map,
    filename = output_file,
    # width = 17.4,      # 单栏宽度（单位：cm）
    # height = 12,
    # units = "cm",
    dpi = 300,
    compression = "lzw",
    bg = "white"
  )
})

cat("All maps saved to:", normalizePath(output_dir))

# 生成各年度地图对象列表
map_list <- map(years, ~{
  plot_cluster_map(
    cluster_data = st_clean_en,
    map_data = map_data,
    target_year = .x
  ) + 
    tm_layout(legend.show = F)  # 隐藏单个图例
})

map_list[[8]]


# 修改第20张图的创建代码（重点调整legend部分）
legend_map <- tm_shape(map_data) + 
  tm_polygons() +
  tm_layout(bg.color = "transparent") +
  tm_compass(#position = c("left", "top"), 
             position = c(0.2,0.8),
             size = 1.5) +
  tm_scale_bar(
    # position = c("right", "top"),
    position = c(0.2,0.6),
               breaks = c(0,400,800),
               text.size = 0.7) +
  tm_add_legend(
    type = "fill",
    labels = c("Primary Cluster", "Secondary Cluster", "Non-cluster"),
    col = c("#d73027", "#fc8d59", "white"),
    title = "Cluster Category",
    border.lwd = 0.5
  ) +
  tm_layout(
    legend.only = TRUE,
    legend.position = c(0.2,0.3),
    legend.width = 1,#1.5,  # 增加图例宽度
    legend.title.size = 0.8,#1.5,
    legend.text.size = 0.8#1.1#,  # 适当减小字号
    # legend.text.width = 1.2  # 允许标签换行
  )
legend_map
# 将说明图加入列表
combined_list <- c(map_list,list(legend_map))

# 组合拼接（5列4行布局）
final_map <- tmap_arrange(
  combined_list,
  ncol = 5, 
  nrow = 4,
  outer.margins = 0.03#,
  # asp = 0.6  # 调整长宽比
)
final_map


# 高质量输出


tmap_save(
  final_map,
  filename = "主图/时空扫描地图/Combined_Map.png",
  width = 45,     # 宽度 45cm（期刊通常接受 10-50cm）
  height = 35,    # 高度 35cm
  units = "cm",   # 单位厘米（期刊推荐）
  dpi = 600       # 分辨率 600 DPI（期刊要求 ≥ 300 DPI）
)
# 
# tmap_save(final_map,
#           filename = "主图/时空扫描地图/Combined_Map.tiff",
#           width = 45, 
#           height = 35, 
#           units = "cm",
#           dpi = 600,
#           compression = "lzw")

# 4.1 附图：逐年感染状态统计图--------
# 数据预处理和绘图参数设置
status_counts <- read_excel("./output/表格5.1_感染状态_新发感染区县数量.xlsx")
yearly_incidence_county <- read_excel( "./output/表格3.1.1_发病数和发病率_逐年.xlsx")
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)

# 1. 基础参数设置 
total_counties <- n_distinct(yearly_incidence_county$county_code)

# 2. 计算年度状态分类 
status_data <- yearly_incidence_county %>%
  group_by(county_code) %>%
  arrange(year) %>%
  mutate(
    cum_cases = cumsum(cases),
    prev_cum = lag(cum_cases, default = 0)
  ) %>%
  ungroup() %>%
  mutate(
    status = case_when(
      cases > 0 & prev_cum == 0 ~ "Newly infected",
      prev_cum > 0 ~ "Previously infected since 2005",
      TRUE ~ "No cases reported"
    ),
    status = factor(status, levels = c("Newly infected", "No cases reported", "Previously infected since 2005"))
  )

# 3. 统计各状态区县数 
status_counts <- status_data %>%
  count(year, status, name = "count") %>%
  complete(year, status, fill = list(count = 0))

# 4. 颜色规范（IDoP期刊风格） 
status_colors <- c(
  "Newly infected" = "#D73027",
  "No cases reported" = "#D9D9D9",
  "Previously infected since 2005" = "#FDAE61"
)

# 5. 学术级堆叠柱状图 
ggplot(status_counts, aes(x = factor(year), y = count, fill = status)) +
  geom_col(
    position = "stack",
    width = 0.75,
    color = "white",
    linewidth = 0.2
  ) +
  geom_text(
    aes(label = ifelse(count > 0, count, "")),
    position = position_stack(vjust = 0.5),
    size = 3.5,
    color = "black"
  ) +
  scale_fill_manual(
    values = status_colors,
    breaks = c("Newly infected", "No cases reported", "Previously infected since 2005"),
    guide = guide_legend(reverse = FALSE, title = NULL)
  ) +
  scale_y_continuous(
    expand = expansion(mult = c(0, 0.02)),
    breaks = seq(0, total_counties, by = 10)  # 根据实际调整刻度
  ) +
  labs(
    x = NULL,#"Year",
    y = "Number of Counties"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    axis.title = element_text(face = "bold"),
    axis.text.x = element_text(
      angle = 45,
      vjust = 1,
      hjust = 1,
      color = "black"
    ),
    legend.position = "top",
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")
  )

# 6. 保存高清图 
ggsave("./主图/感染状态柱状堆叠图（正确）.png", width = 8, height = 6, dpi = 600, bg = "white")


# 5 主图5：诊断延误图--------

## 1.5 诊断-延误：发病至诊断时间情况（中位数，四分位间距）有问题
# 读取数据并初步处理
data <- read_excel("./data/processed/cleaned_data9.xlsx") %>%
  filter(year >= 2005) %>% 
  mutate(birth_date = ymd(birth_date), 
         onset_date = ymd(onset_date),
         diag_time = ymd(diag_time))
psych::describe(data$diag_delay)
summary(data$onset_date)
summary(data$diag_time)

data$county_code <- recode(data$county_code, `652201` = '650502', `654203` = '654223')


#下面这些人的发病时间有问题【待讨论】

data %>% filter(onset_date <= as.Date("2003-01-01")) #9人
summary(data$year);summary(data$diag_delay) 
data_debug <- data %>% filter(onset_date >= as.Date("2005-01-01"))
# data_debug <- data %>% filter(diag_delay <=365)
summary(data_debug$diag_delay) #清理后，最大诊断延误时间为3901天
# data <- data %>% dplyr::select(card_id, onset_date:month) 

# 定义一个函数来移除异常值
remove_outliers <- function(df, column) {
  qnt <- quantile(df[[column]], probs=c(.25, .75), na.rm = T)
  iqr <- IQR(df[[column]], na.rm = T)
  df[!(df[[column]] < (qnt[1] - 5*iqr) | df[[column]] > (qnt[2] + 5*iqr)), ] #20250321改
}

# 应用函数到数据框
data_cleaned <- remove_outliers(data, "diag_delay")
summary(data_cleaned$diag_delay)


# Infectious Diseases of Poverty 双图组合可视化 
library(ggplot2)
library(dplyr)
library(cowplot)

### 数据预处理 
data_cleaned <- data_cleaned %>% 
  mutate(
    year = factor(year, levels = 2005:2023)) %>% 
  filter(diag_delay >= 0 )  # 排除异常值

### 可视化主题设置
idp_theme <- function(base_size=10) {
  theme_minimal(base_family = "Arial") %+replace%
    theme(
      plot.title = element_text(face="bold", size=12, hjust=0.5),
      axis.title = element_text(face="bold", size=10),
      axis.text = element_text(color="#333333"),
      panel.grid.major.x = element_line(color="grey90", linewidth=0.3),
      panel.grid.minor.x = element_blank(),
      legend.position = "none"
    )
}

### 图1：全省年度延误分布 ----
p1 <- ggplot(data_cleaned, aes(x=diag_delay, 
                               y=factor(year, levels=rev(levels(year))))) +  # 逆序因子)) +
  geom_boxplot(
    fill = "#F28E2B",
    outlier.size = 1,
    notch = TRUE  # 显示统计显著性凹槽[1](@ref)
    # fill = "#F28E2B",
    # outlier.shape = 21, 
    # outlier.color = "#E15759",
    # width = 0.7
  ) +
  stat_summary(  # 添加中位数标注[3](@ref)
    fun = median, 
    geom = "text", 
    aes(label = round(..x.., 1)),
    size = 3,
    color = "#333333",
    hjust = -0.3
  ) +
  geom_vline(xintercept = 7, linetype="dashed", color="#D62728") +  # WHO阈值
  scale_x_continuous(
    # name = "Diagnosis delay (days)",
    breaks = seq(0, 60, 15),
    limits = c(0, 60)
  ) +
  labs(
    y = "Calendar year",
    title = "Diagnosis Delays by Year"#,
    # subtitle = "Xinjiang Autonomous Region, 2005-2023"
  ) +
  idp_theme()

### 图2：地级市延误分布 ----
data_cleaned %>% names()
data_cleaned$county_code %>% unique()
nameCE_city <- nameCE_city %>% 
  arrange(nameE)
nameCE_city$nameE <- factor(nameCE_city$nameE)

data_cleaned1 <- data_cleaned %>% 
  mutate(city_code = stringr::str_sub(county_code, 1, 4)) %>% 
  left_join(nameCE_city, by = "city_code") 

data_cleaned1 <- data_cleaned1 %>%
  mutate(nameE = factor(nameE, levels=rev(sort(unique(nameE))))) 

# data_cleaned1$city_code %>% unique();nameCE_city$city_code %>% unique()
# names(data_cleaned1)
data_cleaned1$nameE %>% unique()
 
p2 <- ggplot(data_cleaned1, aes(x=diag_delay, y=reorder(nameE, diag_delay))) +
  geom_boxplot(
    fill = "#F28E2B",
    outlier.size = 1,
    notch = TRUE  # 显示统计显著性凹槽[1](@ref)
  ) +
  stat_summary(  # 添加中位数标注[3](@ref)
    fun = median, 
    geom = "text", 
    aes(label = round(..x.., 1)),
    size = 3,
    color = "#333333",
    hjust = -0.3
  ) +
  geom_vline(xintercept = 7, linetype="dashed", color="#D62728") +  # WHO阈值
  scale_x_continuous(
    # name = "Diagnosis delay (days)",
    breaks = seq(0, 60, 15),
    limits = c(0, 60)
  ) +
  labs(
    y = "Prefecture-level City",
    title = "Diagnosis Delays by City",,
    # caption = "Data source: Xinjiang CDC surveillance system | Median values annotated on right"
  ) +
  idp_theme() +
  theme(axis.text.y = element_text(size=8))

### 双图组合排版 ----
combined_plot <- plot_grid(
  p1, p2,
  ncol = 2,
  labels = c("A", "B"),
  label_fontface = "bold",
  align = "v",
  axis = "lr"
)
combined_plot

### 输出图形
ggsave("./主图/诊断延误图.tiff",
       plot = combined_plot,
       width = 18, 
       height = 25,
       units = "cm",
       dpi = 600,
       compression = "lzw")

# 6 主图6：空间自相关图--------
combined_df <- readRDS("./output/LISA结果_历年.rds")

library(tidyverse)
library(tmap)
library(sf)
library(here)

# 数据预处理 ----------------------------------------------------------------
# 合并地理数据与分析结果
map_data_prep <- map_data %>%
  mutate(county_code = as.character(GB1999)) %>%
  left_join(combined_df, by = "county_code") %>%
  mutate(
    cluster_type = case_when(
      cluster_type == "Low-High" ~ "Cold Spot",
      cluster_type == "High-Low" ~ "Hot Spot",
      cluster_type == "Low-Low" ~ "Cold Cluster",
      cluster_type == "High-High" ~ "Hot Cluster",
      TRUE ~ "Non-significant"
    ) %>% factor(levels = c("Hot Spot", "Hot Cluster", "Cold Spot", "Cold Cluster", "Non-significant")),
    # 筛选显著结果 (p < 0.05)
    cluster_type = if_else(Pr_Ii < 0.05, as.character(cluster_type), "Non-significant")
  )

# 定义地图绘制函数 ----------------------------------------------------------
plot_lisa_map <- function(target_year) {
  # 过滤当年数据
  annual_data <- map_data_prep %>%
    filter(year == target_year)
  
  # 期刊配色方案 (ColorBrewer 8-class RdYlBu)
  cluster_colors <- c(
    "Hot Spot" = "#d73027",    # 深红色
    "Hot Cluster" = "#fc8d59", # 橙色
    "Cold Spot" = "#4575b4",   # 深蓝色
    "Cold Cluster" = "#91bfdb",# 浅蓝色
    "Non-significant" = "gray90"
  )
  
  # 创建地图
  tm_shape(annual_data) +
    tm_polygons(
      col = "cluster_type",
      palette = cluster_colors,
      title = "Cluster Type",
      border.col = "gray80",
      lwd = 0.2,
      alpha = 0.8
    ) +
    tm_layout(
      main.title = paste("Spatial Autocorrelation:", target_year),
      main.title.size = 0.9,
      main.title.position = "center",
      # legend.show = FALSE,  # 隐藏单个图例
      # compass.show = FALSE, # 隐藏指北针
      # scale.show = FALSE    # 隐藏比例尺
      fontfamily = "Arial"

    )
}

# 创建独立图例页 ----------------------------------------------------------
create_legend_page <- function() {
  tm_shape(map_data_prep) + 
    tm_polygons() +
    tm_layout(bg.color = "transparent") +
    tm_compass(
      position = c(0.12, 0.85),
      size = 2.5,
      color.dark = "black"
    ) +
    tm_scale_bar(
      position = c(0.12, 0.78),
      breaks = c(0, 200, 400),
      text.size = 0.8
    ) +
    tm_add_legend(
      type = "fill",
      labels = c("Hot Spot", "Hot Cluster", "Cold Spot", "Cold Cluster", "Non-significant"),
      col = c("#d73027", "#fc8d59", "#4575b4", "#91bfdb", "gray90"),
      title = "Cluster Category",
      border.lwd = 0.3
    ) +
    tm_layout(
      legend.only = TRUE,
      legend.position = c(0.15, 0.5),
      legend.title.size = 0.9,
      legend.text.size = 0.7,
      legend.width = 0.25
    )
}

# 主程序 --------------------------------------------------------------------
# 生成各年度地图
years <- 2005:2023
map_list <- map(years, plot_lisa_map)

# 添加图例页
legend_page <- create_legend_page()
combined_list <- c(map_list, list(legend_page))

# 组合输出 (5列4行布局)
final_map <- tmap_arrange(
  combined_list,
  ncol = 5,
  nrow = 4,
  outer.margins = c(0.02, 0, 0.02, 0),
  asp = 0.6
)

# 高质量输出
tmap_save(
  final_map,
  filename = here("Figures", "Combined_LISA_Maps.tiff"),
  width = 45,       # 单位：厘米
  height = 35,
  units = "cm",
  dpi = 600,
  compression = "lzw"
)

